/**
 * @file   Mesh.h
 * @author ych <ych@Ubuntu18-04>
 * @date   Sun Jun 27 11:22:10 2021
 * 
 * @brief Base class for all types, including irregular and regular.
 * 
 * 
 */

#ifndef __YCH_MESH__
#define __YCH_MESH__

#include "Point.h"
#include <string>
#include <fstream>
#include <iostream>
#include <cstdlib>
#include <vector>


#define TEMPLATE template <size_t DIM, typename T>

TEMPLATE class Mesh;

TEMPLATE class IrregularMesh;

TEMPLATE
Point<DIM, T> mid_point(const Point<DIM, T> &p0, const Point<DIM, T> &p1);

TEMPLATE
class Mesh
{
public:
    virtual size_t get_n_points() const = 0;
    virtual size_t get_n_edges() const = 0;
    virtual size_t get_n_grids() const = 0;
    virtual size_t get_n_dofs() const = 0;
    virtual const Point<DIM, T> &get_point(size_t _idx) const = 0;
    virtual const Geometry &get_edge(size_t _idx) const = 0;
    virtual const Geometry &get_grid(size_t _idx) const = 0;
    virtual std::vector<int> get_all_boundary() {};
    virtual Point<2, double> &get_dof(size_t _idx) = 0;

    // 判断某编号的自由度是否是边界，目前只有结构网格用，非结构网格可以用，等人写
    virtual bool IsBoundary(size_t idx){};

    // 给结构网格的接口，非结构网格不用
    virtual const Point<2, double> &get_lbc() const{};
    virtual const Point<2, double> &get_ruc() const{};
    virtual size_t get_nx() const{};
    virtual size_t get_ny() const{};
    virtual size_t CortoIdx(int i, int j){};
};

class Easymesh : public Mesh<2, double>
{
private:
    std::valarray<Point<2, double>> nodes;
    std::valarray<Geometry> edges;
    std::valarray<Geometry> grids;
    std::valarray<Point<2, double>> dofs;

public:
    void readData(const std::string &filename);
    Easymesh() : Mesh(){};
    Easymesh(const std::string &filename);
    size_t get_n_points() const;
    size_t get_n_edges() const;
    size_t get_n_grids() const;
    size_t get_n_dofs() const;
    const Point<2, double> &get_point(size_t _idx) const;
    const Point<2, double> &get_dof(size_t _idx) const;
    const std::valarray<Point<2, double>> &get_dofs() const;
    const Geometry &get_edge(size_t _idx) const;
    const Geometry &get_grid(size_t _idx) const;
    Point<2, double> &get_point(size_t _idx);
    Point<2, double> &get_dof(size_t _idx);
    std::valarray<Point<2, double>> &get_dofs();
    Geometry &get_edge(size_t _idx);
    Geometry &get_grid(size_t _idx);
    void build_dofs();
};

const Point<2, double> &Easymesh::get_point(size_t _idx) const
{
    return nodes[_idx];
};

Point<2, double> &Easymesh::get_point(size_t _idx)
{
    return nodes[_idx];
};

const Point<2, double> &Easymesh::get_dof(size_t _idx) const
{
    return dofs[_idx];
};

Point<2, double> &Easymesh::get_dof(size_t _idx)
{
    return dofs[_idx];
};

const std::valarray<Point<2, double>> &Easymesh::get_dofs() const
{
    return dofs;
};

std::valarray<Point<2, double>> &Easymesh::get_dofs()
{
    return dofs;
};

const Geometry &Easymesh::get_edge(size_t _idx) const
{
    return edges[_idx];
};

Geometry &Easymesh::get_edge(size_t _idx)
{
    return edges[_idx];
};

const Geometry &Easymesh::get_grid(size_t _idx) const
{
    return grids[_idx];
};

Geometry &Easymesh::get_grid(size_t _idx)
{
    return grids[_idx];
};

size_t Easymesh::get_n_points() const
{
    return nodes.size();
};

size_t Easymesh::get_n_dofs() const
{
    return dofs.size();
};

size_t Easymesh::get_n_edges() const
{
    return edges.size();
};

size_t Easymesh::get_n_grids() const
{
    return grids.size();
};

/// Most contents copy from Robert Li's AFEPack.
void Easymesh::readData(const std::string &filename)
{
    std::cout << "Reading easymesh data file ..." << std::endl;
    std::cout << "\treading node data ..." << std::flush;
    std::ifstream is((filename + ".n").c_str());
    size_t n_node, n_side, n_element;
    char text[64];

    is >> n_node >> n_element >> n_side;
    is.getline(text, 64);
    nodes.resize(n_node);
    for (size_t i = 0; i < n_node; i++)
    {
	size_t ivtx;
	Geometry::bm_t bm;
	is >> ivtx >> nodes[i][0] >> nodes[i][1] >> bm;
	nodes[i].set_index(ivtx);
	nodes[i].set_boundary_mark(bm);
	nodes[i].set_vertex(0, ivtx);
	nodes[i].set_boundary(0, ivtx);
	nodes[i].set_dof(0, ivtx);
    }
    is.close();

    dofs = nodes;

    std::cout << " OK!" << std::endl;
    std::cout << "\treading side data ..." << std::flush;
    is.open((filename + ".s").c_str());
    size_t n_edge;
    is >> n_edge;
    if (n_edge != n_side)
    {
	std::cerr << "in side file: side number error." << std::endl;
	exit(-1);
    }
    edges.resize(n_side);
    for (size_t i = 0; i < n_side; i++)
    {
	size_t iedg, a, b;
	int ea, eb;
	Geometry::bm_t bm;
	edges[i].set_n_vertex(2);
	edges[i].set_n_boundary(2);
	edges[i].set_n_neighbor(2);
	edges[i].set_n_dofs(2);
	is >> iedg >> a >> b >> ea >> eb >> bm;
	edges[i].set_index(iedg);
	edges[i].set_boundary_mark(bm);
	edges[i].set_vertex(0, a);
	edges[i].set_vertex(1, b);
	edges[i].set_boundary(0, a);
	edges[i].set_boundary(1, b);
	edges[i].set_neighbor(0, ea);
	edges[i].set_neighbor(1, eb);
	edges[i].set_dof(0, a);
	edges[i].set_dof(1, b);
    }
    is.close();
    std::cout << " OK!" << std::endl;

    std::cout << "\treading element data ..." << std::flush;
    is.open((filename + ".e").c_str());
    size_t n_ele;
    is >> n_ele;
    if (n_ele != n_element)
    {
	std::cerr << "in element file: element number error." << std::endl;
	exit(-1);
    }
    size_t n_n;
    is >> n_n;
    if (n_n != n_node)
    {
	std::cerr << "in element file: node number error." << std::endl;
	exit(-1);
    }
    size_t n_s;
    is >> n_s;
    if (n_s != n_side)
    {
	std::cerr << "in element file: side number error." << std::endl;
	exit(-1);
    }
    is.getline(text, 64);
    grids.resize(n_element);
    for (size_t i = 0; i < n_element; i++)
    {
	grids[i].set_n_vertex(3);
	grids[i].set_n_boundary(3);
	grids[i].set_n_neighbor(3);
	size_t igrd, v0, v1, v2, b0, b1, b2;
	int n0, n1, n2;
	is >> igrd >> v0 >> v1 >> v2 >> n0 >> n1 >> n2 >> b0 >> b1 >> b2;
	grids[i].set_index(igrd);
	grids[i].set_boundary_mark(0);
	grids[i].set_vertex(0, v0);
	grids[i].set_vertex(1, v1);
	grids[i].set_vertex(2, v2);
	grids[i].set_neighbor(0, n0);
	grids[i].set_neighbor(1, n1);
	grids[i].set_neighbor(2, n2);
	grids[i].set_boundary(0, b0);
	grids[i].set_boundary(1, b1);
	grids[i].set_boundary(2, b2);
	const Point<2, double> &p0 = nodes[v0];
	const Point<2, double> &p1 = nodes[v1];
	const Point<2, double> &p2 = nodes[v2];
	if ((p1[0] - p0[0]) * (p2[1] - p0[1]) - (p1[1] - p0[1]) * (p2[0] - p0[0]) < 0)
	{
	    std::cerr << "vertices of element " << igrd
		      << " is not correctly ordered." << std::endl;
	    // j = g.vertex(2);
	    // g.vertex(2) = g.vertex(1);
	    // g.vertex(1) = j;
	    // j = g.boundary(2);
	    // g.boundary(2) = g.boundary(1);
	    // g.boundary(1) = j;
	}
	grids[i].set_n_dofs(3);
	grids[i].set_dof(0, v0);
	grids[i].set_dof(1, v1);
	grids[i].set_dof(2, v2);
    }
    is.close();

    std::cout << " OK!" << std::endl;
};

Easymesh::Easymesh(const std::string &filename) : Mesh()
{
    readData(filename);
};

void Easymesh::build_dofs()
{
    int dof_index = 0;
    int n_ele = get_n_grids();
    int n_total_dofs = get_n_points() + get_n_edges();
    dofs.resize(n_total_dofs);
    for (int i = 0; i < n_ele; i++)
    {
	int n_local_dofs = 6;
	grids[i].set_n_dofs(n_local_dofs);
	int n_vtx = grids[i].get_n_vertex();
	for (int k = 0; k < n_vtx; k++)
	{
	    Point<2, double> &vtx = get_point(grids[i].get_vertex(k));
	    vtx.set_dof(0, -1);
	}
	int n_bnd = grids[i].get_n_boundary();
	for (int k = 0; k < n_bnd; k++)
	{
	    Geometry &bnd = get_edge(grids[i].get_boundary(k));
	    if (bnd.get_n_dofs() == 2)
	    {
		int a = bnd.get_dof(0);
		int b = bnd.get_dof(1);
		bnd.set_n_dofs(3);
		bnd.set_dof(0, a);
		bnd.set_dof(1, b);
		bnd.set_dof(2, -1);
	    }
	}
    }
    for (int i = 0; i < n_ele; i++)
    {
	int n_local_dofs = grids[i].get_n_dofs();
	int n_vtx = grids[i].get_n_vertex();
	for (int k = 0; k < n_vtx; k++)
	{
	    Point<2, double> &vtx = get_point(grids[i].get_vertex(k));
	    int global_dof = vtx.get_dof(0);
	    if (global_dof == -1)
	    {
		vtx.set_dof(0, dof_index);
		grids[i].set_dof(k, dof_index);
		vtx.set_index(dof_index);
		dofs[dof_index] = vtx;
		dof_index++;
	    }
	    else
		grids[i].set_dof(k, global_dof);
	}
	int n_bnd = grids[i].get_n_boundary();
	for (int k = 0; k < n_bnd; k++)
	{
	    Geometry &edge = get_edge(grids[i].get_boundary(k));
	    int global_dof = edge.get_dof(2);
	    if (global_dof == -1)
	    {
		edge.set_dof(2, dof_index);
		grids[i].set_dof(k + n_vtx, dof_index);
		Point<2, double> &vtx0 = get_point(edge.get_dof(0));
		Point<2, double> &vtx1 = get_point(edge.get_dof(1));
		Point<2, double> vtx = mid_point<2, double>(vtx0, vtx1);
		vtx.set_dof(0, dof_index);
		vtx.set_boundary_mark(edge.get_boundary_mark());
		vtx.set_index(dof_index);
		dofs[dof_index] = vtx;
		dof_index++;
	    }
	    else
		grids[i].set_dof(k + n_vtx, global_dof);
	}
    }
};

class RegularMesh : public Mesh<2, double>
{
protected:
    Point<2, double> lbc; //x0y0
    Point<2, double> ruc; //x1y1
    size_t nx;
    size_t ny;
    double hx;
    double hy;
    size_t n_dof = 3;

public:
    Point<2, double> node;
    Geometry edge;
    Geometry grid;
    Point<2, double> dof;

    RegularMesh() : Mesh()
	{
	    lbc[0] = lbc[1] = 0;
	    ruc[0] = ruc[1] = 1;
	    nx = ny = 2;
	    hx = hy = 0.5;
	};
    RegularMesh(const Point<2, double> &_lbc, const Point<2, double> &_ruc, int _nx, int _ny);
    const Point<2, double> &get_lbc() const { return lbc; };
    const Point<2, double> &get_ruc() const { return ruc; };
    void set_lbc(const Point<2, double> &_lbc) { lbc = _lbc; }
    void set_ruc(const Point<2, double> &_ruc) { ruc = _ruc; }
    size_t get_nx() const { return nx; };
    size_t get_ny() const { return ny; };
    double get_hx() const { return hx; };
    double get_hy() const { return hy; };
    void set_nx(size_t _nx) { nx = _nx; };
    void set_ny(size_t _ny) { ny = _ny; };
    void set_hx(double _hx) { hx = _hx; };
    void set_hy(double _hy) { hy = _hy; };
    size_t get_n_points() const;
    size_t get_n_edges() const;
    size_t get_n_grids() const;
    size_t get_n_dofs() const;
    Point<2, double> &get_point(size_t _idx);
    virtual Point<2, double> &get_dof(size_t _idx);
    Geometry &get_edge(size_t _idx);
    virtual Geometry &get_grid(size_t _idx);
    void _get_grid(size_t _idx);
    std::vector<int> get_all_boundary();
    bool IsBoundary(size_t idx){
        std::vector<int> _BndIndex = this->get_all_boundary();
        return !(std::find(_BndIndex.begin(),_BndIndex.end(),idx)==_BndIndex.end());
    };
    virtual size_t CortoIdx(int i, int j){return (nx + 1) * i + j;};

    //缺少相关实现，等一个有缘人来填空
    virtual const Point<2, double> &get_point(size_t _idx) const{};
    virtual const Geometry &get_edge(size_t _idx) const{};
    virtual const Geometry &get_grid(size_t _idx) const{};
};

RegularMesh::RegularMesh(const Point<2, double> &_lbc, const Point<2, double> &_ruc, int _nx, int _ny) : Mesh(),
													 lbc(_lbc), ruc(_ruc), nx(_nx), ny(_ny)
{
    hx = (ruc[0] - lbc[0]) / nx;
    hy = (ruc[1] - lbc[1]) / ny;
    grid.set_n_dofs(3);
    grid.set_n_neighbor(3);
}

Point<2, double> &RegularMesh::get_point(size_t _idx)
{
    size_t ix = (_idx) % (nx + 1);
    size_t iy = (_idx) / (nx + 1);
    node[0] = ix * hx;
    node[1] = iy * hy;
    node.set_index(_idx);
    node.set_n_dofs(1);
    node.set_n_vertex(1);
    node.set_vertex(0, _idx);
    node.set_boundary_mark((ix == 0 || iy == 0 || ix == nx || iy == ny));
    return node;
};

Point<2, double> &RegularMesh::get_dof(size_t _idx)
{
    return get_point(_idx);
};

Geometry &RegularMesh::get_edge(size_t _idx)
{
    size_t ix = (_idx) % (3 * nx + 1);
    size_t iy = (_idx) / (3 * nx + 1);
    edge.set_index(_idx);
    if (ix < nx)
    {
	edge.set_neighbor(0, ((iy == 0) ? -1 : (iy - 1) * (2 * nx) + 2 * ix));
	edge.set_neighbor(1, ((iy == ny) ? -1 : iy * 2 * nx + 2 * ix + 1));
	edge.set_dof(0, iy * (nx + 1) + ix);
	edge.set_dof(1, iy * (nx + 1) + ix + 1);
	edge.set_vertex(0, iy * (nx + 1) + ix);
	edge.set_vertex(1, iy * (nx + 1) + ix + 1);
    }
    else
    {
	ix -= nx;
	if (ix % 2 == 0)
	{
	    ix /= 2;
	    edge.set_neighbor(0, ((ix == 0) ? -1 : iy * (2 * nx) + 2 * ix - 1));
	    edge.set_neighbor(1, ((ix > nx) ? -1 : iy * (2 * nx) + 2 * ix));
	    edge.set_dof(0, iy * (nx + 1) + ix);
	    edge.set_dof(1, (iy + 1) * (nx + 1) + ix);
	    edge.set_vertex(0, iy * (nx + 1) + ix);
	    edge.set_vertex(1, (iy + 1) * (nx + 1) + ix);
	}
	else
	{
	    ix /= 2;
	    edge.set_neighbor(0, iy * (2 * nx) + 2 * ix);
	    edge.set_neighbor(1, iy * (2 * nx) + 2 * ix + 1);
	    edge.set_dof(0, iy * (nx + 1) + ix + 1);
	    edge.set_dof(1, (iy + 1) * (nx + 1) + ix);
	    edge.set_vertex(0, iy * (nx + 1) + ix + 1);
	    edge.set_vertex(1, (iy + 1) * (nx + 1) + ix);
	}
    }
    edge.set_boundary_mark((ix == 0 || iy == 0 || ix == nx || iy == ny));
    return edge;
};

Geometry &RegularMesh::get_grid(size_t _idx)
{
    _get_grid(_idx);
    return grid;
}

void RegularMesh::_get_grid(size_t _idx)
{
    size_t dof[3];
    size_t en[3];
    if (_idx % 2 == 0)
    {
	size_t ix = (_idx / 2) % nx;
	size_t iy = (_idx / 2) / nx;

	dof[0] = iy * (nx + 1) + ix;
	dof[1] = iy * (nx + 1) + ix + 1;
	dof[2] = (iy + 1) * (nx + 1) + ix;

	en[0] = _idx + 1;
	en[1] = (iy == 0) ? -1 : _idx - 1;
	en[2] = (ix == 0) ? -1 : _idx - 2 * nx + 1;

	grid.set_dof(0, dof[0]);
	grid.set_dof(1, dof[1]);
	grid.set_dof(2, dof[2]);

	grid.set_neighbor(0, en[0]);
	grid.set_neighbor(1, en[1]);
	grid.set_neighbor(2, en[2]);
    }
    else
    {
	size_t ix = ((_idx - 1) / 2) % nx;
	size_t iy = ((_idx - 1) / 2) / nx;

	dof[0] = iy * (nx + 1) + ix + 1;
	dof[1] = (iy + 1) * (nx + 1) + ix + 1;
	dof[2] = (iy + 1) * (nx + 1) + ix;

	en[0] = (iy == ny) ? -1 : _idx + 2 * nx - 1;
	en[1] = _idx - 1;
	en[2] = (ix == nx) ? -1 : _idx + 1;

	grid.set_dof(0, dof[0]);
	grid.set_dof(1, dof[1]);
	grid.set_dof(2, dof[2]);

	grid.set_neighbor(0, en[0]);
	grid.set_neighbor(1, en[1]);
	grid.set_neighbor(2, en[2]);
    }
};

size_t RegularMesh::get_n_points() const
{
    return (nx + 1) * (ny + 1);
};

size_t RegularMesh::get_n_dofs() const
{
    return (nx + 1) * (ny + 1);
};

size_t RegularMesh::get_n_edges() const
{
    return (nx + 1) * ny + (ny + 1) * nx + nx * ny;
};

size_t RegularMesh::get_n_grids() const
{
    return nx * ny * 2;
};

std::vector<int> RegularMesh::get_all_boundary()
{
    std::vector<int> _BndIndex;
    for(int i = 0;i <= nx;i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(get_n_dofs()- i - 1);
    }
    for(int j = 1;j <= ny - 1;j++)
    {
        _BndIndex.push_back(j*(nx + 1));
        _BndIndex.push_back((j + 1)*(ny + 1) - 1);
    }
    return _BndIndex;
};

class RegularMeshP2 : public RegularMesh
{
public:
    RegularMeshP2() : RegularMesh(){};
    RegularMeshP2(const Point<2, double> &_lbc, const Point<2, double> &_ruc, int _nx, int _ny) : RegularMesh(_lbc, _ruc, _nx, _ny){};
    Point<2, double> &get_dof(size_t _idx);
    Geometry &get_grid(size_t _idx);
    size_t get_n_dofs() const;
    std::vector<int> get_all_boundary();
    size_t CortoIdx(int i, int j){return (2 * nx + 1) * i + j;};

};

Point<2, double> &RegularMeshP2::get_dof(size_t _idx)
{
    size_t ix = (_idx) % (2 * nx + 1);
    size_t iy = (_idx) / (2 * nx + 1);
    node[0] = ix * hx;
    node[1] = iy * hy;
    node.set_index(_idx);
    node.set_n_dofs(1);
    node.set_n_vertex(1);
    node.set_vertex(0, _idx);
    node.set_boundary_mark((ix == 0 || iy == 0 || ix == nx || iy == ny));
    return node;
}

Geometry &RegularMeshP2::get_grid(size_t _idx)
{
    size_t dof[6];
    _get_grid(_idx);
    if (_idx % 2 == 0)
    {
	size_t ix = (_idx / 2) % nx;
	size_t iy = (_idx / 2) / nx;

	dof[0] = 2 * iy * (2 * nx + 1) + 2 * ix;
	dof[1] = 2 * iy * (2 * nx + 1) + 2 * (ix + 1);
	dof[2] = 2 * (iy + 1) * (2 * nx + 1) + 2 * ix;
	dof[3] = 2 * (iy + 0.5) * (2 * nx + 1) + 2 * (ix + 0.5);
	dof[4] = 2 * (iy + 0.5) * (2 * nx + 1) + 2 * ix;
	dof[5] = 2 * iy * (2 * nx + 1) + 2 * (ix + 0.5);
    }
    else
    {
	size_t ix = ((_idx - 1) / 2) % nx;
	size_t iy = ((_idx - 1) / 2) / nx;

	dof[0] = 2 * iy * (2 * nx + 1) + 2 * (ix + 1);
	dof[1] = 2 * (iy + 1) * (2 * nx + 1) + 2 * (ix + 1);
	dof[2] = 2 * (iy + 1) * (2 * nx + 1) + 2 * ix;
	dof[3] = 2 * (iy + 1) * (2 * nx + 1) + 2 * (ix + 0.5);
	dof[4] = 2 * (iy + 0.5) * (2 * nx + 1) + 2 * (ix + 0.5);
	dof[5] = 2 * (iy + 0.5) * (2 * nx + 1) + 2 * (ix + 1);
    }
    for(int i=0;i<6;i++)
	grid.set_dof(i,dof[i]);
    return grid;
}

size_t RegularMeshP2::get_n_dofs() const
{
    return (2 * nx + 1) * (2 * ny + 1);
}

std::vector<int> RegularMeshP2::get_all_boundary()
{
    std::vector<int> _BndIndex;
    for(int i = 0;i <= 2 * nx;i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(get_n_dofs()- i - 1);
    }
    for(int j = 1;j <= 2 * ny - 1;j++)
    {
        _BndIndex.push_back(j* (2  * nx + 1));
        _BndIndex.push_back((j + 1)*(2 * nx + 1) - 1);
    }
    return _BndIndex;
}

template <typename T>
class RectangleMesh : public Mesh<2, T>
{
private:
    std::valarray<Point<2, T>> nodes;
    std::valarray<Geometry> edges;
    std::valarray<Geometry> grids;
    std::valarray<Point<2, T>> dofs;

public:
    RectangleMesh() : Mesh<2, T>(){};
    RectangleMesh(const Point<2, T> &_lbc, const Point<2, T> &_ruc, int _nx, int _ny);
    virtual size_t get_n_points() const;
    virtual size_t get_n_edges() const;
    virtual size_t get_n_grids() const;
    virtual size_t get_n_dofs() const;
    virtual const Point<2, T> &get_point(size_t _idx) const;
    virtual const Point<2, T> &get_dof(size_t _idx) const;
    virtual const std::valarray<Point<2, T>> &get_dofs() const;
    virtual const Geometry &get_edge(size_t _idx) const;
    virtual const Geometry &get_grid(size_t _idx) const;
};

template <typename T>
const Point<2, T> &RectangleMesh<T>::get_point(size_t _idx) const
{
    return nodes[_idx];
};

template <typename T>
const Point<2, T> &RectangleMesh<T>::get_dof(size_t _idx) const
{
    return dofs[_idx];
};

template <typename T>
const std::valarray<Point<2, T>> &RectangleMesh<T>::get_dofs() const
{
    return dofs;
};

template <typename T>
const Geometry &RectangleMesh<T>::get_edge(size_t _idx) const
{
    return edges[_idx];
};

template <typename T>
const Geometry &RectangleMesh<T>::get_grid(size_t _idx) const
{
    return grids[_idx];
};

template <typename T>
size_t RectangleMesh<T>::get_n_points() const
{
    return nodes.size();
};

template <typename T>
size_t RectangleMesh<T>::get_n_dofs() const
{
    return dofs.size();
};

template <typename T>
size_t RectangleMesh<T>::get_n_edges() const
{
    return edges.size();
};

template <typename T>
size_t RectangleMesh<T>::get_n_grids() const
{
    return grids.size();
};

template <typename T>
RectangleMesh<T>::RectangleMesh(const Point<2, T> &_lbc, const Point<2, T> &_ruc, int _nx, int _ny)
{
    T hx = (_ruc[0] - _lbc[0]) / _nx;
    T hy = (_ruc[1] - _lbc[1]) / _ny;
    nodes.resize((_nx + 1) * (_ny + 1));
    edges.resize(_nx * (_ny + 1) + (_nx + 1) * _ny + _nx * _ny);
    grids.resize(_nx * _ny * 2);

    for (int ix = 0; ix < _nx + 1; ix++)
	for (int iy = 0; iy < _ny + 1; iy++)
	{
	    int node_idx = iy * (_nx + 1) + ix;
	    nodes[node_idx].set_index(node_idx);
	    /// debug
	    if (ix == 0 || ix == _nx || iy == 0 || iy == _ny)
		nodes[node_idx].set_boundary_mark(1);
	    /// end
	    nodes[node_idx].set_vertex(0, node_idx);
	    nodes[node_idx].set_boundary(0, node_idx);
	    nodes[node_idx].set_dof(0, node_idx);
	    nodes[node_idx][0] = ix * hx;
	    nodes[node_idx][1] = iy * hy;
	}
    for (int ix = 0; ix < _nx; ix++)
	for (int iy = 0; iy < _ny; iy++)
	{
	    int even_ele_idx = (iy * _nx + ix) * 2;

	    grids[even_ele_idx].set_n_vertex(3);
	    grids[even_ele_idx].set_n_boundary(3);
	    grids[even_ele_idx].set_n_neighbor(3);

	    grids[even_ele_idx].set_index(even_ele_idx);
	    int v0 = iy * (_nx + 1) + ix;
	    int v1 = iy * (_nx + 1) + ix + 1;
	    int v2 = (iy + 1) * (_nx + 1) + ix;
	    grids[even_ele_idx].set_vertex(0, v0);
	    grids[even_ele_idx].set_vertex(1, v1);
	    grids[even_ele_idx].set_vertex(2, v2);
	    grids[even_ele_idx].set_n_dofs(3);
	    grids[even_ele_idx].set_dof(0, v0);
	    grids[even_ele_idx].set_dof(1, v1);
	    grids[even_ele_idx].set_dof(2, v2);

	    int odd_ele_idx = (iy * _nx + ix) * 2 + 1;

	    grids[odd_ele_idx].set_n_vertex(3);
	    grids[odd_ele_idx].set_n_boundary(3);
	    grids[odd_ele_idx].set_n_neighbor(3);

	    grids[odd_ele_idx].set_index(odd_ele_idx);
	    v0 = (iy + 1) * (_nx + 1) + ix + 1;
	    v1 = (iy + 1) * (_nx + 1) + ix;
	    v2 = iy * (_nx + 1) + ix + 1;
	    grids[odd_ele_idx].set_vertex(0, v0);
	    grids[odd_ele_idx].set_vertex(1, v1);
	    grids[odd_ele_idx].set_vertex(2, v2);
	    grids[odd_ele_idx].set_n_dofs(3);
	    grids[odd_ele_idx].set_dof(0, v0);
	    grids[odd_ele_idx].set_dof(1, v1);
	    grids[odd_ele_idx].set_dof(2, v2);
	}
    dofs = nodes;
};

TEMPLATE
Point<DIM, T> mid_point(const Point<DIM, T> &p0, const Point<DIM, T> &p1)
{
    Point<DIM, T> mid_pnt;
    for (int i = 0; i < DIM; i++)
	mid_pnt[i] = (p0[i] + p1[i]) * 0.5;
    return mid_pnt;
};

#undef TEMPLATE

#else
// DO NOTHING.
#endif
